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We performed extensive molecular dynamics (MD) simulations, supplemented by Mode Coupling 
Theory (MCT) calculations, for the Square Shoulder (SS) model, a purely repulsive potential where 
the hard-core is complemented by a finite shoulder. For the one-component version of this model, 
MCT predicted [Sperl et al. Phys. Rev. Lett. 104, 145701 (2010)] the presence of diffusion 
anomalies both upon cooling and upon compression and the occurrence of glass-glass transitions. In 
the simulations, we focus on a non-crystallising binary mixture, which, at the investigated shoulder 
width, shows a non-monotonic behaviour of the diffusion upon cooling but not upon isothermal 
compression. In addition, we find the presence of a disconnected glass-glass line in the phase diagram, 
ending in two higher-order singularities. These points generate a logarithmic dependence of the 
density correlators as well as a subdiffusive behaviour of the mean squared displacement, although 
with the interference of the nearby liquid-glass transition. We also perform novel MCT calculations 
using as input the partial structure factors obtained within MD, confirming the simulation results. 
The presence of two hard sphere glasses, differing only in their hard core length, is revealed, showing 
that the simple competition between the two is sufficient for creating a rather complex dynamical 
behaviour. 



I. INTRODUCTION 

In the last decades, a lot of effort has been devoted 
to understand dynamical arrest in soft matter systems. 
The pioneering investigations of hard-sphere (HS) col- 
loids by Pusey and van Megen [TJ |5] have shown that 
the HS glass transition occurs at a colloidal packing frac- 
tion cf> w 0.58. This transition has been interpreted by 
the ideal Mode Coupling Theory (MCT) [3] for the glass 
transition. Despite suffering of a shift of the actual glass 
transition value, MCT provides a good description of the 
experimental data. The mechanism of arrest is explained 
in terms of the so-called 'cage effect' 0], where par- 
ticles at high densities become trapped by their nearest 
neighbours for an increasingly long time. This mecha- 
nism manifests itself in the form of a two-step decay of 
the density auto-correlation functions approaching the 
liquid-glass transition. 

Subsequent investigations have focused on HS colloids 
in which an additional short-range attraction was added 
with the intent of mimicking the effective interaction (i.e. 
depletion) arising between colloids in suspension with 
non- absorbing polymers 5 . A simple square- well (SW) 
attraction can be used to imitate these systems. MCT 
predictions [6HS] for the dynamics of the SW model at 
high densities revealed an intriguing behaviour. Indeed, 
when the range of the well width A is reduced down to 
a few percent of the particle diameter, a reentrant glass 
line is observed in the temperature-concentration phase 
diagram. This results in two different kind of glasses: 
a first glass (named repulsive glass), which is found at 
high temperature T, is the HS glass driven by the pack- 
ing of particles, while a second glass (named attractive 
glass) is observed at low T, when energetic effects are 



dominant and particles remain caged in their attractive 
wells. In between the two glasses, at intermediate tem- 
peratures, a reentrant liquid region occurs. Therefore 
at the same concentration it is possible to go from one 
glass to the other by lowering T and passing through 
a pocket of liquid states arising from the competition 
between energetic and entropic effects occurring at in- 
termediate T |S]. At even higher densities, a glass-glass 
line is observed, which terminates at an endpoint named 
higher order singularity. Associated to these multiple 
glasses and reentrant melting, MCT predicts the occur- 
rence of anomalous dynamics, which results in a loga- 
rithmic (rather than two-step) decay of the density auto- 
correlation functions approaching the endpoint singular- 
ity, as well as in a subdiffusive behaviour of the particles 
mean-squared-displacement (MSD). Most of these pre- 
dictions have been confirmed by several simulations |10P 
I12j and experiments [S"lll3j. In particular, in simulations, 
it turns out to be useful to draw iso-diffusivity lines [14] 
in the phase diagram. These lines maintain the same 
shape of the MCT glass line at all (sufficiently small) 
values of D, so that they provide a useful reference to es- 
tablish whether a reentrance (and eventually associated 
anomalous dynamics) is present. 

The successful predictions of MCT for the SW system 
paved the way for applications of the theory in a wide 
variety of soft matter systems. In particular, MCT has 
been used to describe the arrested behaviour of several, 
purely repulsive systems. Among these are star polymers 
|15) . i.e. long polymer chains anchored onto a central 
core, where the number of chains (arms) varies the soft- 
ness of the particles, bridging HS colloids (in the limit 
of very large arm number) to polymer chains (when the 
arm number is limited to 2). While one-component star 
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polymer solutions only display a glass driven by pack- 
ing of the stars [16j . binary mixtures of stars of different 
arm numbers and sizes have been shown to display mul- 
tiple glassy states through a combined effort of MCT, 
simulations and experiments |17j . Despite indications 
that anomalous dynamics could be present in these sys- 
tems [IS] , a clear evidence from MCT predictions has not 
been provided. In contrast, a recent theoretical study of 
binary size-asymmetric HS mixtures has reported the 
occurrence of higher-order singularities and a variety of 
different glasses [T3] . 

Without necessarily turning to mixtures, it was re- 
cently realised that one-component systems with distinct 
length scales in the interaction potential (the so-called 
core softened models) are also promising candidates for 
detecting thermodynamic and dynamic anomalies [201 - 
|2"5] , Among these, the simplest model is the square shoul- 
der (SS) model, where the hard-core is complemented by 
an additional repulsive corona. This model has been used 
to describe the behavior of some metallic glasses or 
complex materials like micellar [37] or granular systems 
[2"B] . as well as primitive model of silica [53] and water 

Recent MCT calculations reported the existence of 
multiple glass transitions also for the SS system both un- 
der compression and cooling [30] . A peculiar behaviour 
of the SS model, with no counterpart in other investi- 
gated systems, is the prediction of a disconnected glass- 
glass transition with two endpoint singularities for cer- 
tain values of the shoulder width A. Even though a rich 
phenomenology has been predicted for the SS system, 
numerical simulations aiming to confirm this behaviour 
have not been performed so far. In this work we pro- 
vide an extensive and systematic characterisation of the 
SS model by means of event-driven molecular dynamics 
(MD) simulations in order to describe its dynamical be- 
haviour. We examine the one-component system as well 
as a suitably chosen binary mixture which is considered 
in order to avoid crystallisation at high densities and low 
temperatures, and to probe a sufficiently slow dynam- 
ics. The paper is organised as follows. In Section [IT] 
we describe the simulation methods and provide a sum- 
mary of MCT. Then in Section |III| we report our main 
results in four different subsections: in IIII Al we discuss 
the behaviour of the self-diffusion coefficient calculated 
from the simulations and extract an ideal glass line us- 
ing power-law fits of the data; in |IIIB| we compare with 
existing MCT results and perform new calculations for 
the binary mixture currently under study to closely com- 



pare the theoretical results with the simulations; in HI C 
we then search for the existence of the predicted MCT 
higher order singularities; in |IIID| we report results for 
the non-ergodicity parameters obtained from theory and 
simulations to assess the types of the glasses that the sys- 
tem forms at various packing fractions and temperatures. 
Finally in Section [TV] we discuss our findings and provide 
some conclusions and perspectives. 



II. METHODS: SIMULATIONS AND THEORY 

We study a 50 : 50 mixture of N = 2000 particles of 
species A and B interacting via pairwise SS potential 



r < a,. 



Vij(r) 



(Tij < r < (I + A)er 4 
r > (1 + A)cr 



(1) 



where i,j = A, B, <j aa and obb are the particles diame- 
ters (and (tab = (caa + c r ss)/2), Acr^ = 0.15<7jj are the 
shoulder widths, and uq = I is the shoulder height. The 
mass m of both particles is chosen as unit mass, while 
o bb and Mo are the units of length and energy respec- 
tively. T is measured in units of energy (i.e. /cs=l). 



FIG. 1. Square shoulder potential for a generic additive bi- 
nary mixture of species i,j. Here ct;-, = (l/2)((Ti + <Jj) are the 
hard-cores, Aaij are the shoulder widths and Mo = 1 is the 
shoulder height. 



The size ratio between the two species is oaaI^bb = 
1.2. We also study the simple monodisperse version 
of the SS system with the same width. However, the 
one-component system crystallises before it has actually 
entered a sufficiently slowed-down regime of dynamics, 
similarly to what generally observed for one-component 
glass-formers. The introduction of a small asymmetry in 
the size of the two species favors particles rearrangements 
at state points where the one-component easily crystal- 
lizes. This allows us to investigate states that are several 
orders of magnitude slower than those that are possible 
to explore in the monodisperse case. In this way, we get 
as close as possible to the ideal glass line, which is defined 
as the locus of points in the packing fraction-temperature 
state diagram having diffusivity D — > 0. 

We perform event-driven MD simulations of the sys- 
tem as a function of T and packing fraction, defined 
as <j) = (ir/6)(p A a A + Pbv b ), being p l = NJL 3 , L the 
edge of the cubic simulation box and JV, the number of 
particles for each species. Simulations are performed in 
the canonical and microcanonical ensemble. For the de- 
sired packing fraction an initial configuration is generated 
randomly and the particles velocities are extracted from 
a Maxwell-Boltzmann distribution corresponding to the 
desired T. Then, the system is equilibrated by perform- 
ing MD simulations in the canonical ensemble with ap- 
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propriate rescaling of particle velocities. After the equi- 
libration, for each state point investigated, NVE simula- 
tions are performed for times ranging from t = 10 2 for 
the highly diffusive state points, to t = 10 6 for the most 
viscous state points. For all the simulations t is measured 
in units of OB^m/ito) 1 / 2 . 

A first comparison with MCT results is possible by the 
calculation of iso-diffusivity lines, which typically pre- 
serve the shape of the ideal liquid-glass line. These lines, 
along which the self-diffusion coefficient D is constant, 
are evaluated as follows. We first calculate the mean- 
squared displacement (MSD) (r 2 (t)) of the particles at 
several (<p, T) and extract D from its long-time limit be- 
havior, using Einstein relation, 



D = lim 

t— >oo 



Gf 



(2) 



Then we identify state points with the same D and con- 
nect them by iso-D lines. Repeating this procedure for 
lower and lower values of D, typically covering a few 
orders of magnitude in Z), we can extrapolate the D = 0- 
ideal glass line. To do this, we follow MCT predictions 
for D that should go to zero at the ideal glass transition 
with a power-law dependence. This should apply inde- 
pendently on the chosen path, and hence both along an 
isotherm 



D~\</>-<f> g (T) l 7(T) 
and along an isochorc 



D~\T-T g {4>) 



17(0) 



(3) 



(4) 



Here </> g and T g are, respectively, the critical values of 
the packing fraction and temperature at the ideal glass 
transition, while 7 is a non-universal exponent that is 
also determined by the theory. By performing the power- 
law fits, following Eq. ^ and Q, we can then trace the 
locus of points in the (9, T) phase diagram for which 
D — » 0. This line can be directly compared to the MCT 
glass line, as previously done for other systems [17l EJ 
|3"2] . Indeed, MCT usually overestimates the tendency 
to form a glass, so that the two lines (numerical and 
theoretical) are always shifted by a certain amount in 
both T and <f>. However, the shape of the two lines has 
been found, for all previously investigated systems, to 
be identical: this makes possible to establish an effective 
bilinear mapping between the two curves so that they 
scale on top of each other, as it has been done for the SW 
model [33 . In the presence of singular state points, such 
as the MCT higher order singularities [33], the mapping 
procedure allows to estimate their exact location on the 
numerical phase diagram. Indeed, for the SW system, it 
was shown [3T] that one of such singularities does exist by 
performing ad-hoc simulations near this particular state 
point. In the present work, we aim to carry out a similar, 
detailed investigation for the SS system. 

To clarify some issues that will arise below, we briefly 
summarize the main aspects of MCT here. MCT predicts 



the occurrence of a glass transition starting from a set of 
integro-differential equations for the density correlators 
^g(*) = (Pq*(*)Pq(0))/^(?) a t different wave numbers 
q, where S(q) — (p q * (0)p q (0))/iV is the static structure 



factor and p n (t) = J2"=i exp[iq • r 3 (t)}. The MCT equa- 
tions of motions (for Newtonian Dynamics) [35] read, in 
the one-component case, as 



<i> q (t) + n 2 q <i> q (t) + n 2 q 



dt'm q {t-t')$Jt') = (5) 



where fl 2 = q 2 ksT /mS(q) is a characteristic frequency 
and m q = .Fg is the memory kernel. 

Taking the long-time limit of Eq. ^ one obtains, 



fj{i-h) = r q [fk], 



(6) 



where f q = limt_yoo & q (t) is the so-called non-ergodicity 
parameter. J~ q [fk] is the Mode Coupling functional, 
which is bilinear in f q 



f 



d 3 k 
(2tt) 3 



Vq,k/fc/|q-k|) 



(7) 



where 



^ q .k = S(q)S(k)S(\q - k|)ii[q • kc fc +q- (q- k) C | q _ k ,] 2 

(8) 

and Cfc = 1/[1 — pS(k)] is the direct correlation function. 
At the glass transition MCT predicts that, for t —> 00, 
the correlator does not decay to zero but reaches a finite 
plateau value. 

From Eq. Q it is clear that the only inputs needed 
to solve Eq. (TSp are the number density p and the static 
structure factor S(q) of the system. The latter can be 
obtained by solving the Ornstein-Zernike equation [36] 
through the use of integral equations or it can be evalu- 
ated numerically from simulations. 

Besides the occurrence of a liquid-glass transition, un- 
der specific conditions, a system can display multiple 
glassy states, giving rise to the presence of glass-glass 
transitions in the kinetic phase diagram. These multiple 
glasses occur as bifurcations of the solutions of Eq. ^ 
upon variation of the control parameters. Across a glass- 
glass transition the non-ergodicity parameter jumps dis- 
continuously between two non-zero values. This transi- 
tion is found to terminate at an endpoint, named higher- 
order singularity, beyond which one can go from one 
glassy solution to the other continuously. The higher- 
order singularities can be of type A 3l when the two 
glasses coalesce already inside the glassy region, and of 
type A± when the two glasses merge also with the liquid 
solution right on top of the liquid glass line. The latter 
is a very special point occurring at (</>*, T* , A*), which 
can be identified by finely tuning the value of the control 
parameter A [8] [30] , and in its vicinity the form of the 
decay of <& q {t) is predicted to be unique. 

Solving the full dynamical Eq. (p| close to any point 
on the liquid-glass transition, $„(t) is found to follow a 
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typical two-step decay. A first decay at short times cor- 
responds to the characteristic time that particles employ 
to explore the cages formed by their nearest neighbors. 
A second decay occurs at longer time, and is character- 
ized by the a-relaxation time associated to the structural 
rearrangements necessary for restoring the ergodicity in 
the fluid. In between these two regimes, $ g (t) displays 
a characteristic plateau which is associated to the size 
of the cages in which particles are rattling before finally 
escaping. It is well-accepted that the long-time relax- 
ation of the correlators can be described by a stretched 
exponential 

$,(*) ~ ^exp-e/^ (9) 

where f q , r q give an estimate respectively of the non- 
ergodicity parameter, the a-relaxation time, while f3 q is 
the stretching exponent. 

While no general analytic solution of the MCT equa- 
tions is provided for $ 9 (t), its asymptotic form close to 
the glass is known. On approaching the liquid-glass tran- 
sition, the correlators are described by the Von Schwei- 
dler power-law decay [3J. However, close to a higher or- 
der singularity, *& q (t) shows a peculiar logarithmic de- 
pendence: 

%(t) ~ 11 h& Ht/r q ) + h& \n\t/r q ). (10) 

The parameters fq,h q ,h q are the critical non- 
ergodicity parameter and critical amplitudes of first and 
second order in the expansion in ln(t) [3]. It is found that 

a specific value of the wave vector q* exists, at which h q 2 ^ 
is zero, thus allowing for a pure logarithmic decay of the 
correlator to be observed. Hence, the correlators should 
display a characteristic concave (convex) shape for q < q* 
(q > q*) in a logarithmic time scale. 

For the SW system, higher-order singularities have 
been predicted and observed by numerical simulations 
and experiments. In particular, it was shown [51137) that 
MCT predictions in this case are robust upon the use of 
different closure relations such as mean spherical approx- 
imation (MSA) or Percus-Yevick (PY). Different closures 
only produce a shift of the glass transition lines with re- 
spect to each other. 

For the SS system the situation appears to be more 
complex. Recent theoretical studies EE] have shown 
that for the same value of the shoulder width, the use 
of two different closures, namely PY and Rogers- Young 
(RY), as input to the theory (from now on denoted as 
RY-MCT and PY-MCT respectively), provide qualita- 
tively different results. While the liquid-glass line ob- 
tained within RY-MCT displays two reentrances (and 
hence diffusion minima and maxima) associated both to 
cooling (as in the SW system) and to compression, no 
reentrance is observed using PY-MCT. In the latter case, 
there is also no evidence of a glass-glass transition, while 
RY-MCT predicts two glass-glass lines each terminating 
in a higher order singularity. However, for the SS sys- 
tem at the investigated A, RY is expected to be superior 



to PY in the description of its structural and thermody- 
namic properties. Therefore, one of the aims of this work 
will be also to assess the validity of PY-MCT or RY-MCT 
predictions in order to establish the correct scenario for 
the SS system while approaching the glass transition. 



III. RESULTS 

A. Iso-diffusivity lines and ideal liquid-glass line 
from simulations 

We start by reporting the behavior of D along isother- 
mal and isochoric cuts in the (4>,T) phase diagram in or- 
der to assess the presence of diffusion anomalies, perhaps 
like the ones observed in other core-softened potentials 
[2T1 l22l 139] , conceptually similar to the SS system. We 
remark that all shown data points do not crystallise and 
have reached a diffusive behaviour at long times, a con- 
dition necessary in order to extract D. In the following, 
we report results only for A particles, because the be- 
haviour of the B particles is qualitatively the same due 
to the quasi-one-component nature of the mixture. 

Fig. [2ja) shows the normalized self-diffusion coeffi- 
cient of A particles Da/Dq as function of <j> for several 
isotherms, varying from T = 5.0 (close to the HS regime) 
to T = 0.3 (where the shoulder effect becomes promi- 
nent). The normalization factor Dq = (Jbb\JT /m is in- 
troduced to account for the T dependence of the particles 
average velocities. The behavior of Da with cj) is similar 
along all studied isotherms: from the dilute limit Da de- 
creases monotonically at all T. For not too low T, the 
decrease becomes faster with increasing </> and is compat- 
ible with a power-law decay, as discussed below. Notable 
exceptions are data points for T < 0.35: in these cases, 
a robust power-law dependence is not observed with de- 
creasing T. Indeed, at T = 0.3 the data show a much 
wider range of decay. For a deeper investigation we have 
performed simulations for T < 0.3 and many adjacent <f> 
with mesh 0.05 in the range 0.40 < <fi < 0.53. This has 
allowed us to carefully check that the behavior of Da is 
strictly monotonic within numerical error at the studied 
A value of the SS model. However, despite the absence of 
a diffusivity maximum, we are tempted to speculate that 
at these low T the observed slower decrease of Da seems 
to be an effect of the competition between the two length 
scales in the potential. Indeed, at low enough 4> the sys- 
tem behaves as being composed of effective HS particles 
of diameter a + A, while with increasing <f> the bare hard- 
core at a becomes dominant, thereby providing an inter- 
mediate non-trivial (^-dependence of Da- Such scenario 
does not exclude the presence of a non-monotonic behav- 
ior for different values of A, that will be investigated in 
future studies. 

Next, we investigate the behaviour of Da/Dq with T 
along several different isochores, ranging between <fi = 
0.40 and (f> = 0.59. This is reported in Fig. [2j». At 
lower <j) the system does not reach a glass transition even 
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FIG. 2. Normalized diffusion coefficient Da/ Dq as a func- 
tion of (a) cj> for several isotherms, as reported in the labels. 
At all investigated T data show a monotonic decrease with 
increasing <j>, which clearly indicates the absence of diffusion 
anomalies associated to compression/expansion. A crossing of 
the data at high <f) however signals the presence of a diffusivity 
maximum associated to cooling; (b) T for several isochores, 
as reported in the labels. 



for very low temperatures. Indeed, in this limit, it can 
be considered as an effective hard sphere of diameter a + 
A. Since hard spheres are expected to undergo a glass 
transition at (j>^ s ~ 0.58, we expect for the low-T limit 
that the system becomes glassy for <f>^ s [A/(a + A)] 3 ~ 
0.381, in agreement with our simulations. 

For packing fractions 0.40 < <fi < 0.56 the diffusion 
coefficient decreases monotonically with T. For tj> > 0.57, 
Da/ Do becomes non-monotonic: from the high T limit, 
it initially increases and then decreases, giving rise to the 
presence of a diffusivity (local) maximum at intermediate 
T. This is also visible from the crossing of the high-T 
data (at large <j>) in Fig. [2ja) . The anomalous behavior 
of Da upon cooling is similar to that observed for the SW 
system at high enough packing fractions when the width 
of the well is of few percent of the particle size [TT] . 

Compiling all data from Figs. J2fa),(b) we are able to 
trace isodiffusivity lines in the phase diagram to be com- 
pared with the MCT glass lines. The monotonic (non- 
monotonic) behaviour of Da is reflected in the absence 



(presence) of a reentrance in the iso-D^ lines. 

Fig.^a), shows the iso-diffusivity curves for three fixed 
values of normalised diffusion coefficients: Da/Dq = 
1.0 x 10 _3 ,1.0 x 10~ 4 and 1.1 x 10~ 5 . Each curve is 
obtained by extrapolating from Fig. [2ja) and Fig. [2^b) a 
set of i states (<^i,Tj) having the same value of Da/D - 
As discussed above, we do not observe a reentrant be- 
haviour along <f> even if we consider very low values of 
Da/D . On the other hand, for Da/Dq ~ 10~ 5 a reen- 
trance is observed along T, as highlighted in the inset of 
Fig. |3ja). This is in agreement with the presence of a 
diffusivity maximum at high enough <j> (Fig. [2jb)). 

From this analysis we conclude that the liquid-glass 
line shows only a reentrance along T. While the glass 
at high T corresponds to a HS glass, the low-T glass 
could have a different nature and its properties will be 
elucidated in the following. 
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FIG. 3. (a) Isodiffusivity lines for Da/Dq = 1.0 x 1CT 3 , 1.0 x 
10~ 4 and 1.1 x 10 -5 , as well as the extrapolated arrest (Da = 
0) lines from the fits Da ~ \4> — <^g(T)| 7 ' T ' along isotherms 
and Da ~ \T — Tgffi^' along isochores. The data display a 
reentrance in T (inset), while no reentrance in <j> is observed; 
(b) Power law fits along isotherms (left) and isochores (right). 



As said above, the iso-D^ lines are precursors of the 
ideal liquid-glass line, where Da — > 0. We can extrapo- 
late this line where the simulated system should undergo 
dynamical arrest by performing power-law fits of the dif- 
fusivity both along isochores and along isotherms, follow- 
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TABLE I. Extrapolated values of y(T), <f> g , 7O) and T g ob- 
tained from fitting data of Fig. [3] (a) and (b) with MCT pre- 
dictions of Eqs. (3]4l for the diffusion coefficient Da- Error 
bars of the fit parameters typically amount to a few percent 
for the values of 4>g an< A T g , while the 7 exponents can vary 
systematically over different fit intervals, so they should be 
taken with caution. 



ing the MCT predictions in Eq.Q and Eq.Q. The fits 
are shown in Fig. pm). In this way, we can extract both 



the transition values 4> g and T g where the system arrests 
and the associated power-law exponents 7. These values 
are reported in Table [IJ Note that no fits are performed 
for T < 0.35 due to the fact that the data do not show a 
clear power-law behaviour. Also for low values of <f> some 
deviations from power-law dependence are observed at 
low T. 

The resulting ideal glass line is also shown in Fig. [3] 
As expected the two branches, extrapolated by the dif- 
ferent paths, merge continuously in the high-0, low-T 
region of the phase diagram and confirm the shape of 
the isodiffusivity lines. We note that the power-law ex- 
ponents obtained along each isotherm are consistent with 
previous estimates for HS or SW systems. On the other 
hand, for the fits along isochores the 7 exponents are sys- 
tematically lower, at times going below the lowest limit 
predicted by MCT j3] . However, the values of 7 obtained 
from the fits should be taken with caution due to the sig- 
nificant variation of results upon change of the chosen fit 
interval and relative distance to the transition. Nonethe- 
less, the values of <fi g and T g extracted in the same way 
show only little changes (of the order of a few percent), 
and hence they are robust. 

We note that power-law fits along isochores could not 
be performed in most systems with isotropic potentials, 
where the exploration of the low-T region is preempted 
by intervening phase separation. For systems with direc- 
tional interactions where phase separation is suppressed 
by using a limited valence [401 141) , at low T bonding is 
the dominant mechanism of arrest so that the dynamics 
is dominated by an Arrhenius (strong) behaviour [4"2] . 
Here, however, we do not find evidence of an Arrhe- 
nius dependence even at very low T in the investigated 
window of densities, suggesting that the system remains 
power-law (fragile). Indeed, in the SS system tempera- 



ture does not induce bonding, but rather has an effect on 
the excluded volume of particles by changing the effective 
diameter. In this sense, the arrest at low T remains of 
the same kind of the HS glass, so that a similar behaviour 
(fragile) is then expected throughout the phase diagram. 
However, it is then legitimate to ask, given that the na- 
ture of the glass transition remains the same, whether a 
simple competition between two length scales is capable 
to generate higher-order singularities as those predicted 
by MCT and related glass-glass transitions. 



B. Comparison with old and new MCT results: 
role of the input structure factors and mapping to 
simulations 

We now compare the MD simulation results for the 
ideal glass line with MCT predictions. While a mismatch 
of 4> g and T g values is expected for the theoretical and nu- 
merical glass lines, the two should share the same shape, 
as previously observed for a variety of glass-forming sys- 
tems [TH [TBI HI 123 EH- However, when referring to 
the RY-MCT calculations for A = 0.15 [3D], it is im- 
mediate to notice that while the numerical curve shows 
only one reentrance in T, the RY-MCT results display 
two of them, as shown in the inset of Fig. |4j No reen- 
trance is conversely observed for PY-MCTJ33J. Under 
this situation, we cannot perform a consistent mapping 
as previously done for SW systems [3_T], because the dif- 
ference in the shape of the liquid-glass lines cannot be 
taken into account by a simple rescaling procedure. This 
matter thus deserves further investigation. 

For understanding the difference between theory and 
MD results, at first we investigate the reliability of the 
different closures employed for producing the input S(q) 
entering in MCT. In Fig. [| S(q) evaluated within RY 
and PY closures is shown together with that calculated 
directly from simulations of the monodisperse SS system 
for a representative state point. As expected, PY pro- 
vides a rather poor estimate of S q , since the height of its 
peaks and its amplitudes do not agree with those of the 
S q evaluated from MD, while RY reproduces reasonably 
well the simulation results, as previously found for other 
repulsive potentials [UJ [331 [35]. The good agreement 
between RY and MD S(q) is found for all studied state 
points. However, this comparison is limited to the region 
of the phase diagram where the monodisperse system 
does not crystallize. The quality of the input structure 
factors is reflected in the better agreement of RY-MCT 
with the simulation iso-diffusivity lines. It is therefore 
natural from now on, to refer to RY-MCT results as the 
relevant theoretical predictions for the system. 

However, although the shape of the RY-MCT liquid- 
glass line is more similar to the MD results, with at least 
one reentrance recovered, there is still a discrepancy be- 
tween RY-MCT and MD simulations. In fact, both from 
the analysis of D /D and from the iso-diffusivity lines we 
could not detect the presence of a second reentrance (i.e. 
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a diffusion anomaly) along <fi (at fixed low-T). 




FIG. 4. Static structure factors for a monodisperse SS system 
at T = 0.5, cj> — 0.45 calculated by MD simulations as well as 
solving the Ornstein-Zernike equation within Rogers- Young 
(RY) and Percus-Yevick (PY) closures. Inset: MCT results 
for the liquid-glass and glass-glass lines using PY and RY. 

Given this situation, we performed additional MCT 
calculations explicitly incorporating the binary nature of 
the system under investigation. To avoid to rely on a 
certain closure, we have used as inputs to the theory the 
partial structure factors evaluated from MD simulations 
Sfj IM (q). Hence, we have solved the generalized version 
of long-time MCT equations (Eq. [6]) for a binary mixture 
[4"6] on a discretised grid of 1000 wave vectors up to a 
cut-off value of qasB = 65. This value is sufficient for 
the critical non-ergodicity parameters along the liquid- 
glass line to decay to zero. In this way, we determine the 
liquid-glass and (if any) the glass-glass transition, and 
associated non-ergodicity parameters. 

The resulting liquid-glass line is reported in Fig. [5] to- 
gether with the arrest line extrapolated from the fits of 
Da- Despite the expected shift in the control parameters, 
it now appears that the new MCT results for the mixture 
are in full qualitative agreement with the simulation line, 
since the reentrance in <j) is no longer present. We can 
now operate a bilinear transformation, as previously done 
for the SW system [3TJ E3] > to superimpose the MCT re- 
sults onto the glass line obtained from simulations. The 
parameters are chosen via a best fit procedure, giving as 
a result 



<f> -> 1.10460 + 0.0038 
T^ 0.9052T- 0.0111 



(11) 



and the mapped glass lines are shown in Fig. [5] 

MCT calculations predict a 'disconnected' glass-glass 
line, a scenario that was also present in the one- 
component RY-MCT, albeit for lower values of A [30] . 
Luckily, this glass-glass line lies just inside, but very close 
to the liquid-glass line so that signatures of the two A3 
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S S1M -MCT glass-glass 

♦ Tg 

• % 
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FIG. 5. MCT results for the binary mixture under study us- 
ing the static structure factors calculated from simulations as 
input, labeled as Sq IM — MCT liquid glass (filled squares) 
and glass-glass (open squares). Arrest curve drawn from <f> g 
(filled circles) and T g (filled diamonds) obtained from power- 
law fits of Da as in Fig. [3] Mapped MCT lines onto the arrest 
curve: liquid-glass (filled triangles) and glass-glass (open tri- 
angles). Stars are the two predicted higher order singularities 
A3 and Af . 



endpoints are in principle detectable from simulations. 
Through our mapping, we can now estimate the location 
of the two singularities, that will be referred from now on 
as A\ and A$ , indicating respectively the one at lower 
and higher </>. We find A% = {<j> ~ 0.53, T ~ 0.26) and 
- 0.60, T~ 0.49). 
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C. Searching for higher-order singularities 

In this section we investigate the presence of the 
higher-order singularities predicted by MCT in the nu- 
merical phase diagram. To this aim we concentrate on 
distinct paths in the phase diagram that allow us to ap- 
proach closely the two A3 points. We recall however, 
that both points are buried within the glass region, hence 
they are not directly accessible in equilibrium; moreover 
the behaviour of the observables that we examine are in- 
fluenced also by the presence of the nearby liquid-glass 
transition. In the following we will only concentrate on 
species A, but we stress that the qualitative behaviour is 
identical for type B particles. 

We start by discussing the presence of A\ : we examine 
the dynamical behaviour of the system along the isochorc 
4> = 0.525 with decreasing T. We recall that while the 
endpoint should be found at T ~ 0.26, the system be- 



comes glassy according to MCT for T < 0.28 at this 
volume fraction. 

Fig. [fj] shows the MSD for A particles (r AA ) along 
this path. We observe that upon decreasing T the sys- 
tem shows a peculiar slowing down. Indeed, a char- 
acteristic subdiffusive behaviour at intermediate times 
(0.1 < tD < 10) is observed for T < 0.4. Hence, we ob- 
serve a sort of three-step behaviour of the MSD: after the 
ballistic transient, subdiffusion takes place for roughly 
two decades in time, where (r 2 ) ~ t a with a ~ 0.5. At 
long times the typical pattern of glass-forming systems 
takes place: a plateau later followed by long-time diffu- 
sion. Indeed, when T is very low, the system is approach- 
ing the liquid-glass transition, that manifests itself in the 
MSD as the emergence of the plateau. This is observed 
for T < 0.3 and occurs at tD$ ~ 10 2 . The plateau height 
is found to be ~ 0.04ct^ a . Its square root, which pro- 
vides a measure of the cage or localisation length Iq of 
the glass, turns out to be ~ 0.2a aa, roughly twice the 
typical HS cage length (Iq S ~ O.ler). Indeed, the pack- 
ing fraction is significantly smaller than that of the HS 
ideal glass, due to the effect of the shoulder. While the 
long-time behaviour, and associated plateau, is controlled 
by the liquid-glass transition, the additional intermediate 
behaviour, which indicates the presence of a sub-diffusive 
regime, can be associated to the presence of the higher 
order singularity. 

The presence of subdiffusivity is a hint of a closeby 
higher order singularity, but in order to provide a more 
convincing proof of its existence, we now look at the 
behaviour of the density auto-correlations functions. A 
distinctive feature is the presence of a pure logarithmic 
regime for a certain wave- vector q* , where the second- 
order term of the asymptotic expansion in Eq. [10] van- 
ishes. Below and above q* the data should display a 
typical concave-to-convex transition. To visualise this 
behaviour one should be close enough to the A3 point, 
but far enough from the liquid-glass transition in order 
to avoid that the final two-step decay covers most of the 
time-window and preempts the observation of the log- 
arithmic behaviour. We identify at this 4> the optimal 
temperature obeying these requisites as T = 0.375, for 
which we show the collective normalised density auto- 
correlations functions <J>^ (i) as a function of wave- 
vector in Fig. [7|[a). Indeed, at this T we are able to 
identify q*UAA ~ 7 where the decay of the correlators 
is purely logarithmic. Across q* , the concave-to-convex 
transition in the shape of with time is found. In 

Fig. [7]Jb), the T-dependence of the correlators at fixed 
q = q* is shown. It is clear that with further decreasing T 
the system approaches the liquid-glass transition, so that 
the signal of the logarithmic decay gets lost. Indeed, the 
range where the logarithmic dependence (dashed line) is 
valid shrinks upon reducing T. The evidence reported so 
far points to the existence of a higher order singularity in 
the vicinity of the explored path. While we cannot probe 
its exact location, since the system falls out of equilib- 
rium before this can be accessed, it seems to be located 



within, but not too far inside, the glassy region, compat- 
ibly with MCT predictions of a disconnected glass-glass 
line. 




FIG. 6. MSD for a particles as a function of scaled time tDo 
for (f> — 0.525 as a function of T, indicated in the labels. The 
vertical dotted lines indicate as guides to the eye the regime 
of subdiffusive behaviour, which is highlighted by the dashed 
line (oc i 5 ). 
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FIG. 7. The density autocorrelation functions <& q (t) for 
4> — 0.525 as a function of time for (a) several wave 
vectors at T = 0.375. From top to bottom, qoAA = 
1.88, 2.81, 4.68, 5.63, 7.5, 10.32, 13.12, 17.82, 28.13. A concave- 
convex shape transition is observed around q*OAA ~ 7.0 
where the decay of &^ A (t) is almost purely logarithmic; (b) 
several T at fixed wave vector q — q* . From left to right, 
temperatures are T = 1.0, 0.6, 0.5, 0.4, 0.375, 0.35, 0.3, 0.29. 



Next we investigate the presence of the second sin- 
gularity A3 . To approach it, we monitor the isotherm 
T = 0.5 and check the dynamical behaviour with increas- 
ing <j>. We recall that the endpoint should be located at 
<t> ~ 0.6, while the system should become glassy around 
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4> ~ 0.593, as estimated by the power law fits of Da- 
We repeat the same analysis as done for A% and we find 
very similar results (not shown). Indeed, also in this case, 
we observe subdiffusive regime and logarithmic dynamics 
centred around a similar value of q* . Also in this case, 
the interference of the liquid-glass line does not allow us 
to probe the anomalous time window for a significant 
amount of time. To shed light on this point, simulations 
need be performed for a system with a larger A where the 
glass-glass line merges with the liquid-glass line so that 
the anomalous dynamics can be approached in equilib- 
rium. This investigation is currently in progress. 



D. Non-ergodicity parameters from MCT and 
simulations 

In this section, we show the behaviour of the non- 
ergodicity parameters calculated within MCT and ex- 
tracted from the fits of the density auto-correlation func- 
tions in the simulations. 

We start by showing in Fig. [8] the 'critical' non- 
ergodicity parameters f AA (MCT) for the A species 
along the liquid-glass and the glass-glass lines, calculated 
within MCT. The state points at which f AA (MCT) are 
evaluated are shown and numbered in the inset of Fig. 
[8j We find that along the liquid-glass line the theoretical 
f AA (MCT) follows two different types of behaviour: for 
<f> < 0.50, f AA (MCT) at first decreases, showing a shift 
of the peaks to larger wave vectors (from state point 1 to 
3), while for <f> > 0.50 at all temperatures — below and 
above the reentrance — f AA {MCT) maintains its peaks 
and only grows around them in a continuous way (from 
state point 3 to 5). We note that the f AA (MCT) of the 
lowest and highest T (states 1 and 5) corresponds to two 
glasses made, respectively, by effective hard spheres of 
diameter a + A and hard spheres of diameter a. It fol- 
lows that the two f AA (MCT) can be superimposed on 
top of each other by simply readjusting the diameters. 
In between these two limits, the system experiences the 
competition between the two length scales, which results 
in a non-trivial and non-monotonic behaviour of the non- 
ergodicity parameters. The crossover between the two 
regimes arises in a region that cannot be associated with 
any of the two higher order singularities. 

Furthermore, we also show in Fig. [8] the evolution 
of f AA (MCT)'$ for the 'disconnected' glass, occurring 
upon crossing the glass-glass line. In this case the non- 
ergodicity parameters are much larger and longer-ranged 
in q, indicating more tightly-caged glasses, in analogy 
with previous observations of repulsive glasses for star 
polymer mixtures[18). In the present case, the cage size 
is approximately one half of that of the first glass, as it 
can be estimated by the g-range of the non-ergodicity pa- 
rameters. We observe a non-monotonic behaviour when 
moving along the glass-glass line from low to high density, 
similarly to what observed for the liquid-glass line. 

We now extract the non-ergodicity parameters by fit- 
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q a AA 

FIG. 8. Critical non-ergodicity parameters for the A species 
calculated within MCT along the liquid-glass (curves labeled 
from 1 to 5) and along the glass-glass (curves labeled from 6 
to 8) lines. The corresponding state points and their position 
on the MCT lines are reported in the inset: a non-monotonic 
behaviour with increasing (j> is observed for both sets of data. 



ting the density auto-correlation functions close to the 
arrest line via the stretched exponentials of Eq. ([9| . For 
simplicity, again we focus only on species A, since the 
results are formally identical for both species. 

We show in Fig. p\ the behaviour of f AA along the 
liquid-glass line (for the state points labeled in the upper 
inset). We find a behaviour that is remarkably similar to 
the MCT predictions. Indeed, again we find that the non- 
ergodicity parameter shows a non-monotonic behaviour 
(at fixed q) with increasing <p. While at low <f> and low 
T, f AA is compatible with that of a HS with effective 
diameter a + A, it decreases at first with a shift of the 
peaks at larger q for intermediate densities. Then at a 
certain point, that we can roughly estimate as <j) ~ 0.50, 
the peaks do not move further in q while f AA starts to 
increase in magnitude. This continues until it behaves as 
the f q typical for HS of diameter a. This interpretation is 
reinforced by the fact that we can perfectly scale the low- 
4> fq A 011 top of the high </)-one by simply scaling it for 
the effective diameter (see lower inset). In addition, the 
behaviour of the stretching exponents obtained from the 
fits seems to indicate an increase of (3 AA as the system 
moves along the liquid-glass line with increasing (f>. Hence 
glasses found at lower T are considerably more stretched 
{fi ;$ 0-4) than those found at higher (f>. 

Finally we could not estimate numerically the non- 
ergodicity parameters close the disconnected glass-glass 
transition, due to the fact that this is buried inside the 
glass phase. Hence, we do not have a numerical analogue 
of curves (6-8) in Fig. [8j 
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FIG. 9. (a): Non-ergodicity parameters obtained from sim- 
ulations, fitting the density auto-correlation functions with 
stretched exponentials, for the state points reported in the 
upper inset. The behaviour along the liquid-glass line is strik- 
ingly similar to that of MCT predictions, reported in Fig. [5] 
Lower inset: f AA for low <j> (and low T) (state point 1) is 
identical to that of high <j> (state point 6) upon a rescaling by 
the effective diameter a + A; (b) Stretching exponents j3 AA 
obtained from the fits as a function of wave vector for the 
same state points considered in (a). 



IV. DISCUSSION AND CONCLUSIONS 

In this work we have reported an extensive investi- 
gation by simulations and theory of the SS model. De- 
spite the simplicity of the model, its dynamical behaviour 
close to the glass transition had largely remained unex- 
plored, while its thermodynamic phase diagram showing 
the presence of many crystal phases has been addressed in 
a number of publications . The recent predictions 

by MCT of a unique behaviour of the glass transition for 
the SS model [30] is the motivation for the present work. 
Given the amount of effort requested to elucidate the var- 
ious aspects of both liquid-glass transition and glass-glass 
lines, we have focused on a single value of the width of 



the SS model, referring to future studies to address the 
dependence on A of the results presented here. 

Our investigation initially focused on searching for the 
so-called diffusion anomalies, i.e. a local maximum in 
the self-diffusion coefficient, both upon compressing the 
system and upon cooling it, as predicted by MCT. We 
performed the simulations down to a T-range at the limit 
of today computational capabilities and we did not detect 
the presence of an anomaly upon compression/expansion 
in contrast to theoretical results, as well as with other 
studies of core-softened potentials [HJ [531 (Ml ISO] • How- 
ever, none of these previous studies involved a potential 
with sharply-defined length scales such as the SS, so that 
this might be the reason why a different behaviour is ob- 
served. Also, it remains to be established whether such 
different models share the intriguing glassy properties of 
the one-component SS model, as detailed MCT studies 
have not been performed. To rule out in the SS model 
the presence of a diffusivity maximum in density we will 
need to explore other values of A in the future. Our pre- 
liminary investigations in this direction have also been 
negative in this respect, but we believe that we need to 
perform additional simulations for significantly slow state 
points in order to resolve this controversial aspect. 

The role of the SS width is also crucial within MCT, 
because it controls the position of the glass-glass line and 
its endpoints with respect to the liquid-glass line. The 
present case was predicted to show two glass-glass lines 
merging from the liquid-glass one by RY-MCT for the 
monodisperse SS j3U]. However, we have repeated the 
MCT calculations for the binary mixture investigated by 
MD simulations, which avoided the onset of crystallisa- 
tion, and found a different topology. Namely, a discon- 
nected glass-glass transition lying all inside the glassy 
region with two endpoints, i.e. a scenario that RY-MCT 
attributed to lower values of A. Such a shift in the con- 
trol parameters is expected when comparing with MCT, 
although for the SW model there was a remarkable good 
agreement between theory and simulations at the same A 
[3Tj . Also no reentrance in (j> was detected within the bi- 
nary MCT calculations, in agreement with the simulation 
results. The disconnected glass-glass topology seems to 
be confirmed by simulations, which have shown the pres- 
ence of a regime of logarithmic decay of the correlators 
and of subdiffusive increase of the MSD, in agreement 
with the predictions close to a higher order singularity. 
Also the presence of the glass-glass line and its special 
endpoints, called A% and , is compatible with the 
presented results. However, to really confirm the pre- 
dictions, it will be desirable in the future to approach 
these points in equilibrium, as previously done for the 
SW model, by a fine tuning of A. Again, we plan to 
undertake these investigations in the near future. 

Despite the absence of clear results about the glass- 
glass transition, we have found evidence that the non- 
ergodicity parameters along the liquid-glass line vary in 
a non-monotonic fashion, in very good agreement with 
the theory. We have shown that the glass found at 
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low/intermediate <j> (e.g. <fi — 0.40) and low T is identi- 
cal to that found at high (f> (and almost T-independent) 
when a scaling of the effective diameter is performed to 
take into account the effect of the shoulder width. Thus 
the system clearly displays two identical glasses, both 
of HS type, driven solely by excluded volume, and dif- 
fering between themselves only by a change of length 
scale. However, the interplay between these two glasses 
is highly non-trivial, giving rise to anomalous behaviour 
in the dynamics, although not in the form of a local max- 
imum in the diffusion coefficient, but in the form of clear 
subdiffusivc regime in the MSD and logarithmic dynam- 
ics of the density correlators. This is enhanced at a length 
scale that is compatible with that always associated as 
the main responsible for the hard-sphere transition, i.e. 
the nearest-neighbour length. It is a remarkable find- 
ing that this simple physics is capable of producing such 
non-trivial and unexpected dynamical behaviour. These 
results clearly show that, in order to find higher order 
MCT singularities, it is not needed to rely on two differ- 
ent physical ingredients, such as attraction and repulsion, 
but competing isotropic repulsions are sufficient. 

Finally another aspect that is very peculiar to the 



present study is the fact that at low 4> we can approach 
the glass transition down to very low temperatures, with- 
out intervening crystallisation or phase separation. This 
was previously achieved in other systems with directional 
attraction, where however the dynamics at low T was 
found to be dominated by bonding processes, showing an 
Arrhenius dependence. Here on the contrary, bonding is 
not present and the dynamics seems to always retain the 
fragile character of HS systems. The absence of phase 
separation at low T and <p for the present system thus 
offers the unexplored possibility to investigate the glass 
transition at T — > for <\> — > <f> g [a/(a + A)] 3 for which 
current work is underway. 
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